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Abstract. — It is thought that speciation in phytophagous insects is often due to colonization of novel host plants, because 
radiations of plant and insect lineages are typically asynchronous. Recent phylogenetic comparisons have supported this 
model of diversification for both insect herbivores and specialized pollinators. An exceptional case where contemporaneous 
plant-insect diversification might be expected is the obligate mutualism between fig trees (Ficus species, Moraceae) and 
their pollinating wasps (Agaonidae, Hymenoptera). The ubiquity and ecological significance of this mutualism in tropical 
and subtropical ecosystems has long intrigued biologists, but the systematic challenge posed by >750 interacting species 
pairs has hindered progress toward understanding its evolutionary history. In particular, taxon sampling and analytical 
tools have been insufficient for large-scale cophylogenetic analyses. Here, we sampled nearly 200 interacting pairs of fig and 
wasp species from across the globe. Two supermatrices were assembled: on an average, wasps had sequences from 77% of 6 
genes (5.6 kb), figs had sequences from 60% of 5 genes (5.5 kb), and overall 850 new DNA sequences were generated for this 
study. We also developed a new analytical tool, Jane 2, for event-based phylogenetic reconciliation analysis of very large data 
sets. Separate Bayesian phylogenetic analyses for figs and fig wasps under relaxed molecular clock assumptions indicate 
Cretaceous diversification of crown groups and contemporaneous divergence for nearly half of all fig and pollinator lineages. 
Event-based cophylogenetic analyses further support the codiversification hypothesis. Biogeographic analyses indicate that 
the present-day distribution of fig and pollinator lineages is consistent with a Eurasian origin and subsequent dispersal, 
rather than with Gondwanan vicariance. Overall, our findings indicate that the fig-pollinator mutualism represents an 
extreme case among plant-insect interactions of coordinated dispersal and long-term codiversification. [Biogeography; 
coevolution; cospeciation; host switching; long-branch attraction; phylogeny.] 



Processes affecting the diversification of insects are 
crucial to understanding the origin of biodiversity, 
because most animals are either insect herbivores, 
or natural enemies (predators or parasitoids) of 
these phytophages (Novotny et al. 2002). As primary 
consumers, most insect herbivores are involved in 
antagonistic interactions with plants and, although 
herbivores often exhibit host-specific coevolutionary 
adaptations to plant defenses (Ehrlich and Raven 1964), 
recent empirical studies have suggested that host 
plant lineages are generally older than their associated 
herbivores (Percy et al. 2004; Tilmon 2008; McKenna 
et al. 2009). Such patterns of asynchronous plant-insect 
diversification are consistent with the general paradigm 
that insect speciation results from colonization of novel 



host plants and subsequent reproductive isolation (Percy 
et al. 2004; Tilmon 2008; McKenna et al. 2009; Fordyce 
2010). 

Phytophagous insects are often enemies of plants, 
but some engage in beneficial pollination mutualisms. 
A charismatic example involves the ca. 750 species 
of figs (Ficus, Moraceae) and their pollinating wasps 
(Hymenoptera, Chalcidoidea, Agaonidae) (Fig. 1). 
Agaonids are the only pollen vectors for fig trees and 
agaonid larvae feed exclusively on the flowers of their 
Ficus hosts. Each partner is thus entirely dependent on 
the other for reproduction. Figs are also a major resource 
for frugivores and most animal-dispersed tropical tree 
species interact with vertebrates that also consume 
figs (Howe and Smallwood 1982). The fig-pollinator 
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Figure 1. Classification and worldwide distribution of Ficws. The numbers of species per subgenus is represented as a proportion of total Ficus 
species richness. Breeding systems are indicated as either monoecious (M) or dioecious (D) and modes of pollination are indicated as passive 
(P) or active (A). *Agaon, Alfonsiella, Allotriozoon, Courtella, Elisabethiella, Nigeriella and Paragaon. **Deilagaon and Waterstoniella. 



mutualism is therefore ecologically important in most 
tropical ecosystems (Shanahan et al. 2001). Many fig 
species reproduce irregularly are relatively inaccessible 
in the forest canopy or today are found only in rainforest 
remnants, such that coordinated sampling of Ficus 
and pollinator species for systematic study is difficult. 
These sampling challenges, coupled with the limitations 
of analytical tools for large data sets, have hindered 
progress toward understanding the global evolutionary 
history of the mutualism, despite the fact that many 
details of this intricate symbiosis were described almost 
a century ago (Janzen 1979; Wiebes 1979; Weiblen 2002; 
Cook and Rasplus 2003; Herre et al. 2008). 

Species-specificity in fig pollination appears to be 
extreme compared with most other insect pollination 
mutualisms. Most fig species are pollinated by only 
one or a few wasp species and most wasps are 
associated with just a single fig species (Cook and 
Rasplus 2003; Molbo et al. 2003; Cook and Segar 
2010). Pollinators are specifically attracted to volatile 
compounds emitted by figs (Hossaert-McKey et al. 
1994) and access to the specially modified inflorescences 
is by means of distinctive mandibular appendages 
and detachable antennae (van Noort and Compton 
1996). Pollination is either active (two-thirds of the fig 
species) or passive (one-third, mostly within subgenera 
Pharmacosycea, Ficus, Synoecia, and Urostigma) (Kjellberg 
et al. 2001). Active agaonid wasps collect pollen from 
the anthers of their native figs and store it in thoracic 



pollen pockets (Galil and Eisikowitch 1968; Ramirez 
1978). Once inside a receptive fig, they remove pollen 
from their pockets and deposit it on the flower stigma 
each time they lay an egg (Galil and Eisikowitch 1968; 
Kjellberg et al. 2001). Passively pollinated figs produce 
large quantities of pollen through anther dehiscence and 
wasps are covered with pollen (Galil and Neeman 1977) 
before flying away from their natal figs. 

Closely matching fig and pollinator traits might 
be products of coadaptation (Ramirez 1974; Wiebes 
1979, 1982a; Kjellberg et al. 2001; Weiblen 2004) but, 
regardless, trait-mediated interactions have the potential 
to simultaneously affect the evolution of reproductive 
isolation among pollinator and fig populations; this is 
because fig wasps breed exclusively in pollinated figs. 
This line of reasoning has underpinned the hypothesis 
that cospeciation might account for patterns of fig 
and pollinator diversity. However, this notion runs 
contrary to the paradigm that insect speciation generally 
involves host-switching (Tilmon 2008) and so it remains a 
controversial proposition that requires rigorous testing. 

Under the cospeciation scenario, phylogenies of 
figs and pollinators are expected to show substantial 
congruence. There is some evidence for this pattern 
(Herre et al. 1996; Machado et al. 2005; Ronsted et al. 
2005; Cook and Segar 2010; Cruaud et al. 2011a), but 
recent studies have countered the underlying case for 
cospeciation with evidence of cryptic wasp species and 
relaxed partner specificity. At least 50 fig species are now 
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known to have multiple pollinator species (Michaloud 
et al. 1985, 1996; Rasplus 1996; Kerdelhue et al. 1997; 
Lopez-Vaamonde et al. 2002; Greeff et al. 2003; Molbo 
et al. 2003; Haine et al. 2006; Moe and Weiblen 2010; Chen 
et al. 2012) and as many as 4 different wasp species are 
known to pollinate a single fig species (Machado et al. 
2005; Cook and Segar 2010). Such cases occur in a broad 
taxonomic and geographic spectrum, although cases of 
pollinator species sharing multiple fig species have been 
reported mostly from monoecious figs in the Neotropics 
(Molbo et al. 2003) and the Afrotropics (Erasmus et al. 
2007; Cornille et al. 2012; McLeish and van Noort 2012). 
In any event, evidence of relaxed host specificity and 
some incongruent fig-pollinator phylogenies (Machado 
et al. 2005) suggest that host shifting is a viable 
alternative explanation for fig-pollinator diversification. 

Cospeciation has been hypothesized for the vertically 
transmitted endosymbionts of insects (e.g., Moran 2001; 
Jousselin et al. 2009) but this is not a plausible general 
model for the evolution of plant-insect associations, 
which are horizontally transmitted and not so integrated 
metabolically. Further, if the plant traits that mediate 
insect associations happen to be phylogenetically 
conserved, then host shifting among close relatives 
could also result in topologically congruent phylogenies 
(Percy et al. 2004). In addition, historical biogeography 
has the potential to confound the explanation of such 
patterns if synchronous plant-insect dispersal to new 
environments, followed by geographic isolation, results 
in cospeciation. 

Another useful approach is to investigate patterns 
of temporal congruence (Page and Charleston 1998). 
Divergence time estimates for fig and pollinator clades 
are expected to be approximately equal in the event of 
coradiation, whereas insect lineages are expected to be 
younger than hosts in the case of host shifting (Percy 
et al. 2004; Tilmon 2008; McKenna et al. 2009). 

Previous comparisons of fig and pollinator phylogeny 
have yielded rather different insights on the relative 
importance of host shifting and codiversification 
depending on the taxonomic scope of sampling (Cook 
and Segar 2010). Molecular phylogenetic trees appear 
roughly parallel when based on exemplars of Ficus 
sections and wasp genera (Herre et al. 1996; Jackson 2004; 
Cruaud et al. 2011a), but such deep taxonomic sampling 
is unlikely to detect host shifts among close relatives 
(Machado et al. 2005). On the other hand, regional 
comparisons of particular fig and pollinator clades have 
tended to reject cospeciation in favor of host-switching 
(Machado et al. 2005; Marussich and Machado 2007; 
Jackson et al. 2008; Jousselin et al. 2008), although not 
always (Weiblen and Bush 2002; Silvieus et al. 2008). A 
global test for codiversification therefore requires dense 
sampling of many fig and pollinator lineages across the 
entire geographic range, but a problem of this magnitude 
poses a further methodological challenge. 

Tests of cophylogenetic hypotheses often employ 
tree reconciliation methods that infer evolutionary 
processes such as cospeciation, host shifts, duplications, 
and losses to account for topological incongruence 



between host and associate phylogenies (Page 1994). 
This approach has the power to model the relative 
contributions of different evolutionary processes to a 
given phylogenetic pattern, but biologically realistic 
scenarios become computationally intractable for large 
numbers of taxa (Merkle and Middendorf 2005; Ovadia 
et al. 2011). Genetic algorithms that incorporate dynamic 
programming to efficiently locate and evaluate samples 
from an extremely large universe of event-based 
solutions hold promise in this regard (Conow et al. 
2010). 

Here, we extended the application of a genetic 
algorithm to event-based tree reconciliation analysis for 
cophylogenetic problems involving > 100 taxon pairs and 
applied randomization tests involving null models to 
test the codivergence hypothesis on an unprecedented 
scale. Nearly 200 pairs of interacting fig and fig wasp 
species were sequenced at 5 fig loci (providing up to a 
total of 5.5 kb DNA sequence) and 6 wasp loci (up to a 
total of 5.6 kb). Two supermatrices were assembled. On 
an average, wasps had sequences from 77% of 6 genes, 
figs had sequences from 60% of 5 genes, and overall, 
we generated 850 new DNA sequences for the purpose 
of this study. Maximum likelihood (ML) analyses of fig 
and wasp data sets and Bayesian phylogenetic analyses 
under relaxed molecular clock assumptions enabled 
the comparison of distance, event-based, and temporal 
congruence. Inferences from historical biogeography 
based on our global sample of fig and pollinator clades 
provided additional insight on the relative roles of 
dispersal and vicariance with respect to alternative 
hypotheses of diversification. 



Materials and Methods 
Taxonomic Sampling and DNA Sequencing 

Ficus — We sampled 200 fig species (>l/4 of the 
circa 750 described species) that represent all Ficus 
sections recognized by Berg and Corner (2005) 
(Appendix SI in the Supplementary Material Online, 
doi: 10.5061 /dryad.hr620). Four taxa belonging to the 
tribe Castilleae s.l., Antiaropsis decipiens, Castilla elastica, 
Poulsenia armata, and Sparattosyce dioica, were included 
as outgroups (Datwyler and Weiblen 2004; Ronsted et al. 
2005; Zerega et al. 2005; Clement and Weiblen 2009). Total 
genomic DNA was extracted from 20-30 mg of dried leaf- 
fragments or herbarium material following Ronsted et al. 
(2008). Ficus phylogeny was reconstructed using 5 genes: 
ITS (891 bp), ETS (528 bp), glyceraldehyde 3-phosphate 
dehydrogenase (G3pdh, 769 bp), chloroplast expressed 
glutamine synthetase region (ncpGS, 1630 bp), and 
granule-bound starch synthase (waxy region, 1734 bp). 

Amplification of ITS, ETS, and G3pdh was performed 
following Ronsted et al. (2008). The ncpGS region 
(Emshwiller and Doyle 1999) was amplified using 
Moraceae-specific primers 3F (5/ GTT GTG ATT WAC 
CAT GCT) and 4R (3/ AGA TTC AAA ATC GCC TTC) 
designed for this study. Amplification of ncpGS consisted 
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of 4 min at 94°C followed by 36 cycles of: 1 min 
denaturation (94°C), 1 min annealing (50°C), and 2-min 
extension (72°C). After the last cycle, the temperature 
was kept at 72 °C for a final 5-min extension and then 
lowered to 4°C. The GBSSI or waxy region (Mason- 
Gamer et al. 1999; Clement 2008) was amplified using 
Moraceae-specific primers 3F (5/ GAT CGY GTG TTT 
GTR GAT CAC C) and 10R (3/ GCA ACT GAA TGA 
GAC CAC A). Amplification of waxy consisted of 3 min 
at 94°C followed by 2 cycles of 94°C for 1 min, 58°C for 
1 min, 72°C for 2 min, 2 cycles of 94°C for 1 min, 56°C 
for 1 min, 72°C for 2 min, 2 cycles of 94°C for 1 min, 
54°C for 1 min, 72°C for 2 min, 2 cycles of 94°C for 1 min, 
50°C for 1 min, 72°C for 2 min, and 24 cycles of 94°C 
for 1 min, 48°C for 1 min, 72°C for 2 min. After the 
last cycle, the temperature was kept at 72°C for a final 
20-min extension and then lowered to 4°C. Amplified 
products were purified with the Qiagen PCR purification 
kit (Qiagen Inc.) following the manufacturer's protocols. 
ITS, ETS, G3pdh, and ncpGS were sequenced directly 
from PCR products whereas waxy was cloned using the 
TOPO-TA PCR cloning kit (Invitrogen, Carlsbad, CA). 
Nine clones were screened for inserts, and plasmids were 
isolated from 3 of these using the Qiagen plasmid prep 
kit. Multiple copies of waxy are known in the Rosales 
(Evans et al. 2000) and therefore it was necessary to 
ensure that phylogeny reconstruction was performed 
with orthologous copies. Two copies have been detected 
in Moraceae, GBSSI and GBSS2, which were easily 
distinguished on the basis of size and intron alignment 
(Silvieus et al. 2008; Clement W., unpublished data). 
Analyses were based solely on GBSSI because GBSS2 
was encountered less commonly in figs. 

Cycle sequencing reactions were carried out following 
Ronsted et al. (2008). For sequencing of the ncpGS 
region, internal primers IF (5/ TCW TGW GCT GAA 
AAG CAT), 2F (TTT AAT CTC CAG ACT CSA), and 
5F (5/ TAG TTC ACT CTA AAG GGT) were designed 
for this study in addition to the primers used for 
amplification. Some 50% of the sequences were obtained 
from de novo sequencing for the purpose of this study 
and have been deposited in GenBank (Appendix SI in 
the Supplementary Material Online). Other sequences, 
mostly deposited by coauthors, were obtained from 
existing databases. 

Agaonidae. — Ninety-three percent of the 200 wasps and 
figs used in this study are true associates; i.e., even 
if they were not collected together simultaneously, the 
agaonid species is the pollinator of the fig species. In 
the very few cases where the corresponding agaonid 
was not available in our collection, we used instead 
the pollinator of a closely related species of fig 
(Appendix S2 in the Supplementary Material Online). 
This was always a wasp species that was a close 
congener of the actual pollinator. As the phylogenetic 
position of Agaonidae within the large and complex 
superfamily Chalcidoidea is as yet unknown (Gibson 
et al. 1999; Munro et al. 2011), 4 divergent members 



of the superfamily served as outgroups: Sycophaga 
(Sycophaginae), Ficomila (Eurytomidae), Megastigmus 
(Torymidae), and Trichogramma (Trichogrammatidae). 

All material was collected alive in the field and fixed 
in 95% ethanol. With very few exceptions, Agaonidae 
sequences were obtained from the nondestructive 
extraction of a single wasp specimen (corpse kept 
as voucher). DNA was extracted from a single 
individual that was incubated at 56°C overnight 
(with gentle "shaking" steps by inverting the tubes) 
and using the Qiagen DNeasy kit following the 
manufacturer's protocol. When destructive extraction 
was used, vouchers were selected among specimens 
sampled from the same tree and the same fig after 
careful identification by J.Y.R., Sv.N, and R.U. Vouchers 
are deposited at CBGP, Montferrier-sur-Lez, France. 
To infer phylogenetic relationships between agaonid 
species, we combined 2 nuclear protein-coding genes 
[F2 copy of elongation factor-la (EFla, 516 bp), Wingless 
(Wg, 403 bp)]; 2 mitochondrial protein-coding genes 
[cytochrome c oxidase subunit I (COI, 1536 bp), 
cytochrome b (Cyt b, 749 bp)], and 2 ribosomal genes 
[28S rRNA (D2-D3 and D4-D5 expansion regions, 
1520 bp), 18S rRNA (variable regions V3-5, 787 bp)]. 
Extraction, amplification, and sequencing protocols 
follow Cruaud et al. (2010) for CytB, COI (barcode 
fragment), Wg, 28S, and 18S rRNA, Weiblen (2001) for 
COI [Cl-J-2183 (Jerry)— TL2-N-3014 (Pat) fragment], and 
Cruaud et al. (2011b) for EFla. Both strands for each 
overlapping fragment were assembled using Geneious 
v5.4.2 (Drummond et al. 2007). 

Sixty-seven percent of the sequences were obtained 
from de novo sequencing for the purpose of this study 
and have been deposited in GenBank (Appendix S2 in 
the Supplementary Material Online). Other sequences, 
mostly deposited by coauthors, were obtained from 
public databases. 



Phylogeny Reconstruction 

Protein-coding genes and hypervariable regions were 
aligned using ClustalW 1.81 with the default settings 
(Thompson et al. 1994). Alignments of protein-coding 
genes were translated to amino acids using Mega 4 
(Tamura et al. 2007) to detect frameshift mutations 
and premature stop codons, which may indicate the 
presence of pseudogenes. Alignment of sequences 
encoding rRNA was based on secondary structure 
models (Gillespie et al. 2006), following Cruaud et al. 
(2010). Phylogenetic trees were estimated using both ML 
and Bayesian methods. We selected separate models 
of molecular evolution for different genomic regions 
including mitochondrial genes, rRNA stems, rRNA 
loops + regions of ambiguous alignment, and individual 
nuclear genes using the Akaike information criterion 
implemented in MrAIC.pl 1.4.3 (Nylander 2004). 

For each data set, we performed ML analyses and 
associated bootstrapping (1000 replicates) using the 
MPI-parallelized RAxML 7.0.4 software (Stamatakis 
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2006b). GTRCAT approximation of models (Stamatakis 
2006a) was used for ML boostrapping (1000 replicates). 
Bootstrap percentage (BP) >70% was considered as 
strong support (Felsenstein and Kishino 1993). Bayesian 
analyses (BA) were conducted using a parallel version 
of MrBayes 3.1.1 (Huelsenbeck and Ronquist 2001). 
We assumed across-partition heterogeneity in model 
parameters by considering the parameter m (Nylander 
et al. 2004; Marshall et al. 2006; McGuire et al. 
2007). Parameter values for the model were initiated 
with default uniform priors and branch lengths were 
estimated using default exponential priors. 

To improve mixing of the cold chain and avoid it 
converging on local optima, we used Metropolis-coupled 
Markov chain Monte Carlo (MCMCMC) simulation with 
each run including a cold chain and 3 incrementally 
heated chains. The heating parameter was set to 0.02 in 
order to allow swap frequencies from 20% to 70%. For 
both figs and pollinators, we ran 2 independent runs of 30 
million generations. All the values were sampled every 
3000 generations. For the initial determination of burn- 
in, we examined the plot of overall model likelihood 
against generation number to find the point where the 
likelihood started to fluctuate around a constant value. 
Convergence of the chains was evaluated using the 
online application AWTY (Nylander et al. 2008) and 
the results were based on the pooled samples from 
the stationary phases of the 2 independent runs. Given 
that posterior probabilities (PP) may overestimate clade 
support, for reasons discussed elsewhere (Suzuki et al. 
2002; Cummings et al. 2003; Erixon et al. 2003; Simmons 
et al. 2004), only clades with PP > 0.95 were considered 
strongly supported. All analyses were conducted on 
a 150core Linux Cluster at CBGP, Montferrier-sur-Lez, 
France. 

Test of alternative hypotheses. — To assess whether certain 
alternative relationships among recovered clades could 
be statistically rejected, we performed AU (Shimodaira 
2002) and SH (Shimodaira and Hasegawa 1999) tests in 
the CONSEL package (Shimodaira and Hasegawa 2001). 
The program makermt was used to generate K = 10 sets 
of bootstrap replicates (rl = 0.5, r2 = 0.6, r3 = 0.7, rA = 0.8, 
r5 = 0.9, r6 = l,r7 = 1.1, r8 = 1.2 r9 = 1.3, rlO = 1.4). Each set 
consisted of 100 000 replicates of the row sums (10 times 
the default number of replicates). RAxML was used to 
compute the per-site log likelihoods for all topologies 
tested. To assess the relative support for competing 
phylogenetic hypotheses, we also conducted AU and SH 
tests on recently published data sets (Lopez-Vaamonde 
et al. 2009; Cruaud et al. 2010), which placed Tetrapus as 
sister to all other agaonids with strong support (PP = 0.99 
and BP = 55/59; PP = 1.00, respectively). 

Effects of missing data. — There is debate in the literature 
as to the effect of missing data on the accuracy of 
phylogenetic analyses. Simulation results based on 
limited numbers of characters (Lemmon et al. 2009) 
indicated that nonrandom distributions of missing data 



can result in strong support for nodes that share 
no supporting characters. However, other empirical 
and simulation studies have concluded that taxa with 
extensive missing data can be accurately placed in 
phylogenetic analyses, and that adding characters with 
missing data is generally beneficial, if the overall number 
of characters is large and data are analyzed with 
appropriate methods (see Wiens and Morrill 2011). To 
assess the impact of missing data on our analyses, we 
performed 2 sets of additional analyses. 

First, we built new ("complete species") trees using 
only the more completely sequenced taxa (figs with 
more than 3 genes; wasps with more than 5 genes). 
Then, we used AU and SH tests to compare the full 
(all species and genes) tree, pruned of incompletely 
sequenced taxa, with the matching "complete species" 
tree. Second, we built new ("complete genes") trees by 
removing gene fragments for which <60% of the taxa 
were available. We then used AU and SH tests to test 
if the full tree differed significantly from the "complete 
genes" tree. Taxa were pruned from the combined ML 
tree using the APE package (Paradis et al. 2004) in R 2.14.0 
(http://www.R-project.org). 



Bayesian relative rate tests and long-branch attraction 
artifact. — We tested constancy of evolutionary rates 
among agaonid species using both BEAST 1.5.3 
(Drummond and Rambaut 2007) (coefficient of variation 
statistic and average rate for each branch of the 
chronogram, see "Molecular Dating" section) and a 
Bayesian relative rate (BRR) test (Wilcox et al. 2004). 
For the BBR test, the PP distributions of lengths for 
all branches from the most recent common ancestor 
(MRCA) of the ingroup to each of the terminal taxa 
were based on 1000 randomly chosen post-burn-in 
trees from the BA of the mitochondrial (mtDNA), 
nuclear (nuDNA) and combined data sets, respectively. 
Following Wilcox et al. (2004), we considered rates of 
evolution significantly different between 2 taxa if the 95% 
confidence interval of the PP distribution of the summed 
branch length did not overlap. Branch length estimates 
were compiled using Cadence vl.08b (Wilcox et al. 2004). 

Long-branch attraction (LBA) artifacts can be difficult 
to detect, but methods have been proposed and 
we applied these to our data (Bergsten 2005). For 
computation time reasons, all additional analyses were 
conducted using RAxML only. 

Removing first and third codon positions, which are 
fast evolving, may be a way to reduce LBA. However, 
this can also compromise tree resolution (Kallersjo et al. 
1999; Savolainen et al. 2002; Stefanovie et al. 2004). The 
RY-coding strategy (Woese et al. 1991), by discarding 
fast-evolving transitions and reducing compositional 
bias, constitutes a better approach (Phillips and Penny 
2003; Philippe et al. 2005). We therefore compared the 
topologies obtained with or without RY-coding of: (i) 
the third (nt3) and (ii) first (ntl) and third mtDNA codon 
positions. 
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Long-branch extraction is another approach 
advocated for cases where LBA is suspected (Pol 
and Siddall 2001). Because LBA to the outgroups is 
the most frequent problem, analyses were conducted 
without the outgroups (Bergsten 2005). 

Finally, the different sensitivity of parsimony and ML 
methods can help to detect if LBA is playing a major 
role (Brinkmann et al. 2005). We therefore performed 
parsimony analysis on our data set to detect potential 
shifts in position of agaonid groups. Parsimony analyses 
were conducted with TNT version 1.1 (Goloboff et al. 
2008), using New Technology Search: 1000 replicates 
of random addition sequences, followed by random 
sectorial searches with default options, 100 cycles of 
ratchet and 3 rounds of tree-fusing. All substitutions 
were equally weighted and gaps treated as missing 
data. Robustness of topologies was assessed by bootstrap 
procedures using 1000 replicates. 



Cophylogenetic Analyses 

We tested the congruence between fig and wasp 
phylogenies using both distance and event /topology- 
based methods. The former generate patristic distance 
matrices between species in each phylogeny and 
then test for correlations between the 2 matrices. In 
contrast, event-based methods use evolutionary events 
[cospeciation, duplication, host-shifts, lineage sorting, 
and "failure to diverge" (Page and Charleston 1998; 
Charleston and Perkins 2006)] to map the associate 
phylogeny to the host one. A cost is assigned to each 
event type and we seek to find mappings that minimize 
the total cost. Statistical analyses can be performed by 
comparing the best costs found for the host-parasite data 
set against those of randomized instances. 

We used the distance-based method, ParaFit, 
developed by Legendre et al. (2002) and implemented 
in the program CopyCat (Meier-Kolthoff et al. 2007). 
ParaFit evaluates the global hypothesis of host- 
associate cospeciation with a matrix permutation test of 
codivergence. This test combines 3 types of information: 
the associate phylogeny and the host phylogeny both 
described by their respective matrices of patristic 
distances, and the observed host-associate links. Each 
matrix representing associates and hosts is transformed 
into a matrix of principal coordinates. The association 
is then described by a new matrix, which includes 
both matrices of principal coordinates and the matrix 
of association. Patristic distances were computed from 
fig and wasp ML-phylogenetic trees. Tests of random 
association (null hypothesis) were performed using 
9999 permutations globally across both phylogenetic 
trees. Although the distance-based approach is 
computationally simple, it only yields a measure of 
overall phylogenetic congruence and no information 
on the relative distribution of underlying evolutionary 
events that might have produced the pattern. 

Event-based methods have the advantage of modeling 
evolutionary processes directly, but are computationally 



intensive (Charleston 2009). The problem of finding a 
mapping (event-based reconstruction) of minimum total 
cost has been shown to be computationally intractable 
("NP-complete") (Ovadia et al. 2011). Some existing 
software packages, e.g., TREEMAP 1.0; (Page 1994) 
and 2.02; (Charleston and Page 2002), use exhaustive 
searches, which are prohibitively slow and also permit 
only limited numbers of species. Other programs use 
heuristics (Merkle and Middendorf 2005), which are fast 
but may converge on suboptimal or invalid solutions 
(e.g., ancestral speciation inferred to have occurred 
after speciation of descendants nodes). For this reason, 
our analyses used a genetic algorithm to search a 
sample of the possible solution space with a dynamic 
programming step that efficiently evaluates the cost 
of each such sample. This approach, which finds 
solutions of near-optimal cost, was first implemented 
in the Jane software package (Conow et al. 2010) 
and was validated using a number of existing data 
sets in the literature (Libeskind-Hadas and Charleston 
2009). However, the sheer size of our data sets put 
the analysis far beyond the computational limits of 
the original version of Jane, which also lacks support 
for randomization tests. We therefore substantially 
optimized and improved the existing Jane cophylogeny 
software package, resulting in a new system, Jane 2, 
which is capable of performing event-based analyses 
of very large data sets. Jane 2 and its tutorial are 
freely available for research and educational purposes 
at http:/ /www.cs.hmc.edu/~hadas/jane/index.html. 

We used Jane 2 with the following parameter values: 
the number of "generations" (iterations of the algorithm) 
was set to 40 and the "population" (number of samples 
per generation) was set to 1000. We explored 3 different 
cost models, each with 2 types of randomization test. 
The first model set costs per event as cospeciation = 0 
and all other events = 1. This correponds to the 
TreeMap cost scheme so that a duplication event actually 
contributes two to the total cost because each of the 
2 daughter lineages contributes one duplication event. 
The cost of the best solution was compared with 
the costs found in 100 randomizations in which the 
tip mappings were permuted at random, a method 
advocated by Aldous (2001). The second randomization 
involved 100 randomly generated pollinator trees, of 
the same size as the actual wasp pollinator tree, with 
random tip mappings. The random pollinator trees were 
constructed using the Yule model with beta parameter 
equal to —1. 

In the second model, we used costs of 0 for 
cospeciation, 1 for each duplication, 1 for each host 
switch, and 2 for each loss event. In the third model, 
we set the cospeciation cost at -1 and all other costs to 0, 
where a negative cost maximizes the number of inferred 
cospeciations. For the second and third cost models, we 
used the same 2 randomization tests described for the 
first model. 

All the analyses were performed at Harvey Mudd 
College (Claremont, CA) on a heterogeneous cluster of 
commodity computers comprising a total of 168 cores. 
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On a single commodity computer (e.g., a dual core 
Macintosh), a single fig/ wasp tree required ~3 h of 
computation and thus 100 randomized trials required 
several hours on our cluster. 



Molecular Dating 

We used the uncorrelated log-normal relaxed clock 
method implemented in BEAST 1.5.3 (Drummond and 
Rambaut 2007) and the same modeling strategies as 
for MrBayes and RAxML analyses. We assumed a Yule 
tree prior and we used default priors for all other 
parameters. We used 2 runs of 60 million generations 
with sampling every 6000 generations for figs, and 2 
runs of 240 million generations with sampling every 
24 000 generations for wasps. The 2 separate runs 
were then combined using LogCombiner 1.5.3. We 
ensured convergence using TRACER 1.5 (Drummond 
and Rambaut 2007). Following the removal of 10% 
burn-in, the sampled posterior trees were summarized 
using TreeAnnotator 1.5.3 to generate a maximum clade 
credibility tree and calculate the mean ages, 95% highest 
posterior density intervals (95% HPD) and PR We used 
independent calibration points to estimate divergence 
ages of the main Ficus and agaonid clades. Following 
Ronsted et al. (2005), crown-group Ficus was assigned 
a uniform prior distribution with a minimum age of 
60 Ma based on fossilized achenes (Collinson 1989) 
and a maximum age of 198 Ma based on converging 
molecular estimates for the origin of the angiosperms 
(Bell et al. 2005). Given uncertainties over the age 
of Dominican amber (Grimaldi 1994; Iturralde-Vinent 
and MacPhee 1996, 1999), crown-group Pegoscapus and 
Tetrapus were assigned uniform prior distributions with 
minimum ages of 15 Ma and maximum ages of 60 Ma 
based on Dominican amber fossil (Poinar 1993; Penalver 
et al. 2006). For both fig and wasp phylogenies, nodes 
including taxa endemic to La Reunion were modeled 
with a normal distribution with a mean of 8 Ma 
and SD of 1 Ma based on the proposed age for the 
Mascarene archipelago (McDougall and Chamalaun 
1969; McDougall 1971). 



Ancestral Area Reconstructions and Evolution of Pollination 

Mode 

We inferred the evolution of pollination mode and 
the ancestral areas for figs and their pollinators using 
both ML and parsimony approaches implemented 
in Mesquite 2.73 (Maddison and Maddison 2008). 
Pollination modes and ancestral areas were inferred 
on the ML topologies. For ML optimization, we used 
a stochastic Markov model of evolution (Mkl). The 
Likelihood Decision Threshold was set to 2 log- 
likelihood units. Character data for Ficus and Agaonidae 
were obtained both from the literature (Kjellberg 
et al. 2001; Berg and Corner 2005) and from our 
examination of flowers, pollen pockets, and coxal combs. 



Following Lopez-Vaamonde et al. (2009), current species 
distributions were categorized into 4 character states: 
(0) Afrotropics, (1) Australasia, (2) Neotropics, (3) 
Eurasia. However, because several taxa occur in both 
Eurasian and Australasian regions and a couple of 
taxa occur in both Eurasian and Afrotropical regions, 
and Mesquite requires unique character states, we also 
defined 2 other states: (4) Australasia + Eurasia and 
(5) Afrotropics + Eurasia. We took into account all 
published geographic localities for Ficus and agaonids, 
museum specimens and about 3000 samples of fig wasp 
communities that we collected over the last 15 years. We 
also used the dispersal-extinction-cladogenesis model 
implemented in Lagrange (Ree and Smith 2008), using 
the same raw data and 4 character states. Dispersal rate 
between all areas was set to 1 during the whole period 
considered (data available upon request). 



Results 
DNA Sequence Data 

The completeness of taxa in the combined data 
matrices is different for fig and wasps (Appendices 1-2 in 
the Supplementary Material Online and Supplementary 
Table SI). On an average, wasps have sequences from 
77% of the 6 genes and 67% of the species were sequenced 
for at least 5 gene regions. On an average, figs have 
sequences from 60% of the 5 genes and 70% of the 
species were sequenced for at least 3 regions. Plastid 
regions provide little phylogenetic information within 
Ficus, enforcing the use of more informative single copy 
nuclear regions. These are known to be notoriously 
difficult to amplify from plants in general (Ronsted et al. 
2007) and this was also the case for Ficus in the present 
study. Indeed, ncpGS and waxy matrices only include 
24 and 23% of the taxa, respectively. Models chosen by 
MrAIC for each partition were as follows. Ficus data set: 
GTR + T (ETS, ITS, ncpGS, and waxy), GTR + 1 + Y (G3pdh); 
Agaonid data set: GTR + T (mtDNA), GTR + 1 + r (EFla, 
Wg, rRNA stems), HKY + I + r (rRNA loops). Given 
that a and the proportion of invariable sites can not 
be optimized independently from each other (Gu 1995) 
and following Stamatakis' personal recommendations 
(RAxML manual), we used GTR + r with 4 discrete 
rate categories for all partitions. As RAxML does not 
implement the HKY model, we used GTR instead. 



Wasp Phylogeny 

Our phylogenetic trees (Fig. 2 and Supplementary 
Fig. SI), reconstructed using ML and Bayesian 
approaches provide several new insights into the 
systematics of fig wasps. 

Monophyly of the genera and intergeneric relationships. — 
Fifteen agaonid genera are recovered as monophyletic 
with strong support (Agaon, Alfonsiella, Allotriozoon 
Ceratosolen, Courtella, Deilagaon, Elisabethiella, Eupristina, 
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Kradibia, Nigeriella, Pegoscapus, Pleistodontes, Tetrapus, 
Valisia, and Waterstoniella). In contrast, Platyscapa is 
polyphyletic, and Dolichoris is paraphyletic with respect 
to Blastophaga psenes (the pollinator of F. carica and type 
species of the genus Blastophaga), indicating the need for 
taxonomic rearrangements (Cruaud et al. 2012). 

The relationships among the major clades are unclear 
(Fig. 2). BEAST analysis places Ceratosolen + Kradibia 
(subfamily Kradibiinae) as the sister group to 
the remaining Agaonidae with strong support 
(PPbeasT = 0.98), but this position is not strongly 
supported by MrBayes (PPMrBayes = 0-88) and ML 
analyses (BP = 43). BA place Tetrapus (monogeneric 
subfamily Tetrapusinae) nested within the Agaonidae 
With Strong support (PP M rBayes = 1-00, PPbeast = 1-00), 
although this position is only moderately supported by 
ML analyses (BP = 67). 

Phylogenetic placement of the genus Tetrapus. — By not 
placing Tetrapus as sister to all other agaonids, our 
topology challenges all previous molecular studies by 
ourselves and others (Herre et al. 1996; Machado et al. 
1996, 2001; Lopez-Vaamonde et al. 2009; Cruaud et al. 
2010; Supplementary Table S3). This result deserves 
further examination, so we have conducted additional 
analyses on not only our current data set, but also 
previously published data sets. We provide here a 
summary of the main results (see the Appendix S3 in the 
Supplementary Material Online for further details): 

(1) Both AU and SH tests fail to reject alternative 
topologies in which either Tetrapus, or the clade of 
pollinators associated with subgenus Synoecia and 
subsection Frustescentiae (corresponding to Group 4 
in Cruaud et al. 2010), is constrained to be the sister 
group to all other Agaonidae (Supplementary Table 
S2). Furthermore, AU and SH tests also fail to reject 
alternative positions of Tetrapus using 2 previously 
published data sets (Lopez-Vaamonde et al. 2009; 
Cruaud et al. 2010) that recover Tetrapus as sister to 
all other Agaonidae (Supplementary Table S2). 

(2) Tetrapus is recovered nested within the Agaonidae 
in all the analyses conducted to assess the impact of 
missing data on the accuracy of our phylogeny. AU 
and SH tests showed that phylogenetic trees pruned 
of incompletely sequenced taxa and trees built only 
on gene fragments for which at least 60% of the 
taxa were available, were not significantly different 
from the original trees including all available data 
(Supplementary Fig. S2B,C and Supplementary 
Table S2). Therefore, our analyses show that missing 



data are not responsible for the position of the genus 

Tetrapus. 

(3) Examination of branch lengths (mtDNA, nuDNA, 
and combined trees, Supplementary Figs. SI and 
S3) indicates considerable variation in rates of 
molecular evolution among agaonid lineages This 
result is confirmed by the BRR tests (Supplementary 
Figs. S4 and S5) and the BEAST outputs (95% 
credible interval for the coefficient of variation of 
rates is not abutting against zero for each partition 
and covariance values span zero). Furthermore, a 
long branch leading to Tetrapus, is visible in both 
the nuDNA tree (Supplementary Fig. S3b) and 
ML and Bayesian combined trees (Supplementary 
Fig. SI). BRR tests and branch-specific rates inferred 
by BEAST reveal a lineage-specific increase in 
nucleotide substitution rates on this branch, and 
this is also the case for the branch leading to the 
outgroups (Supplementary Fig. S4). 

(4) RY-coding of first and third mtDNA codon positions 
does not result in significant topological changes, 
but increases support for Tetrapus nested within 
the Agaonidae (Supplementary Fig. S2D,E and 
Supplementary Table S2). 

(5) Unrooted and rooted topologies appeared congruent 
(Supplementary Fig. S2J), showing that rooting does 
not alter the ingroup topology. Furthermore, the 
unrooted topologies from Lopez-Vaamonde et al. 
(2009) and Cruaud et al. (2010) do not show conflicts 
with the topology presented here (Supplementary 
Fig. S7b,c). 

(6) Parsimony analysis of the combined data set recovers 
Tetrapus as sister to the remaining Agaonidae 
(BP = 64) (Supplementary Fig. S6). 

We conclude that neither our study nor previous ones 
have a strong basis for inferring which group is sister 
to all other agaonids. Accordingly, the placement of 
Tetrapus remains unresolved. However, we suggest that 
the repeated recovery of Tetrapus as sister to all other 
agaonids in previous studies may be due to LBA to the 
outgroups and we await further studies. 



Ficus Phylogeny 

The Ficus phylogenetic trees (Fig. 2 and 
Supplementary Fig. S9) are globally congruent with 
previous hypotheses (Herre et al. 1996; Weiblen 2000; 
Jousselin et al. 2003; Ronsted et al. 2005, 2008; Cruaud 
et al. 2011a; Xu et al. 2011) (Supplementary Table S5). 



Figure 2. BEAST chronograms of the evolutionary history of figs and fig wasps. Groups of figs and their associated genera of pollinators 
are represented using the same color. Ficus subgenera and Agaonidae subfamilies according to current classifications are delimited by colored 
rectangles (Pharma. for Ficus subgenus Pharmacosycea). Pie charts at main nodes show the likelihood of different geographic areas of origin as 
inferred by Mesquite (see "Materials and Methods" section). Gray rhombuses show clades of fig species from Continental Asia, whereas gray 
arrows indicate hypothesized southward migration of clades. Squares correspond to node supports: Black square: BP > 70% and PPMrBayes 
or PPbeast > 0.95; white square: BP > 70% or PPMrBayes or PPbeast > 0.95. Details in this figure can be viewed at greater magnification at 
Systematic Biology online. 
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Monophyly of the subgenera and infrageneric relationships. — 
Several moderately to strongly supported clades 
broadly correspond to currently recognized sections or 
subsections based on previous molecular phylogenetic 
studies (see Ronsted et al. 2008) and morphology 
(sections Pharmacosycea, Oreosycea, Americana, 
Galoglychia, Adenosperma s.L, Sycomorus s.L, Sycocarpus, 
Eriosycea, and subsections Malvanthera, Conosycea, 
Urostigma, Ficus, and Frutescentiae). Only 3 of the 6 Ficus 
subgenera currently recognized based on morphology 
(Berg and Corner 2005) are recovered as monophyletic 
with strong support. These are: Sycomorus (BP = 71, 
PPMrBayes = 0.75, PP B east = 1-00); Sycidium (BP = 100, 
PPMrBayes = 1-00, PPbeast = 1-00); and Synoecia (BP = 100, 
PPMrBayes = 1-00, PPbeast = 1-00). Relationships of 
deeper nodes are not strongly supported. The first split 
within Ficus is between section Pharmacosycea (BP = 100, 
^MrBayes = 1-00, PPbeast = 1-00) and the remainder 
of Ficus (BP = 39, PP M rBayes = 0.88, PPbeast = 0.85). 
The next split is between a clade with all members 
of subgenus Urostigma except subsection Urostigma 
(BP = 100, PPMrBayes = 1-00, PP B EAST = 1-00) and a 
clade with members of subsection Urostigma, subgenus 
Sycomorus and all other dioecious figs (BP = 66, 
PPMrBayes = 0.95, PP B EAST = 1-00). 

Exploration of bias in the Ficus phylogenetic trees. — Previous 
molecular studies are similar in recovering section 
Pharmacosycea (pollinated by the genus Tetrapus) as sister 
to the other Ficus species. However, with the exception of 
the BEAST analysis by Xu et al. (2011), this relationship 
is supported by parsimony only (Supplementary Table 
S5). The difference in likelihood scores between our 
best ML tree and the trees from analyses constrained to 
place either subgenus Sycomorus or a clade of subgenera 
{Sycomorus, Sycidium, Ficus, and Synoecia) sister to the 
remaining Ficus were not significant (Supplementary 
Table S2). This confirms that relationships within Ficus 
are unstable along the backbone of the tree and should 
be regarded as uncertain. 

Analyses conducted to assess the impact of missing 
data on the accuracy of our phylogeny resulted in 
topologies that were congruent with the topology 
estimated from the global data set (Supplementary Table 
S2). It is noteworthy that using only Ficus species for 
which at least 3 gene regions were available slightly 
increases node support, but deeper nodes remain 
unresolved (not shown). 



Cophylogenetic Comparisons 

All our analyses rejected the null hypothesis of no 
correlation between fig and wasp phylogenies. Using 
distance-based methods, the global test of cospeciation 
(Parafit) rejected a random association between host 
and pollinator taxa (ParaFitGlobal = 1.37866, P < 0.01). 
Further, 176 of the 200 tests of individual host-associate 
pairs resulted in significant associations between figs 



and their agaonid pollinators (P < 0.01) (Supplementary 
Table S6). 

In event-based analyses, exact results depend on 
the weights assigned to different speciation events. 
Under the classic TreeMap cost-model of zero for 
cospeciation and one for other events (Charleston 
and Page 2002), Jane 2 inferred 198 cospeciation 
events, 204 duplications, 102 host shifts and 61 losses 
between fig and wasp phylogenies, accounting for an 
optimal cost of 367. Whatever the cost model used, the 
number of cospeciation events inferred by Jane 2 was 
always significantly greater than expected by chance 
(Supplementary Fig. S10). 

This topological correlation suggests codiversifi- 
cation, but does not establish a time line, so we next used 
independent relaxed molecular-clock dating techniques 
to test for contemporaneous divergence (Percy et al. 
2004). We found strong temporal congruence between 
both stem and crown mean ages of most partner clades 
and between the ages of inferred cospeciation events 
(Fig. 3). Codiversification test results were not sensitive 
to the order of deep branching in the phylogenies. 



Evolution of Pollination Modes 

Parsimony and likelihood reconstruction on the 
wasp topology both inferred the ancestral pollination 
mode as ambiguous. Parsimony inferred passive and 
active pollination as equiprobable ancestral conditions. 
Similarly, the likelihood difference between the 2 states 
was not significant (proportional likelihoods of 0.53 
and 0.47, respectively) (Supplementary Fig. Sll). Using 
the Ficus topology, parsimony again inferred active 
and passive pollination as equibrobable. However, 
likelihood favors active pollination as the ancestral 
condition (proportional likelihood of 0.91 versus 0.09 for 
passive pollination). Overall, the reconstructions reveal 
that pollination modes are homoplastic with several 
independent shifts between states (passive/ active) along 
both phylogenies (Supplementary Fig. Sll). 



Biogeographic Analyses 

Our dating analyses indicate that the current 
pantropical distribution of the mutualism cannot have 
resulted simply from vicariance following the break- 
up of Gondwanaland. Instead, our ancestral area 
reconstructions suggest that figs and their pollinators 
arose simultaneously in Eurasia (Mesquite proportional 
likelihood = 0.72 for figs and 0.97 for wasps, 
Supplementary Fig. S12) during the Late Cretaceous 
~75 Ma (74.9 Ma for figs and 75.1 Ma for wasps; Fig. 2, 
Supplementary Fig. S13, and Table 1). Mesquite and 
Lagrange results were similar, indicating that fig wasps 
most probably arose in Eurasia. However, Lagrange 
reconstructions for the fig phylogeny were equivocal 
due to a basal polytomy. The Eurasian region was 
proposed as the ancestral area of origin for Ficus in one 
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Figure 3. Temporal evidence for fig and fig wasp codivergence. a) Correlation between stem and crown mean ages of major fig and wasp 
groups (with 95% HPD). b) Temporal congruence of the 198 cospeciation events inferred by Jane 2. 



of the alternative reconstructions that fall within 2 log- 
likelihood units of the optimal scenario (data not shown, 
but available upon request). Although the concordance 
in means crown ages is striking, the PP density around 
the mean estimate is quite wide (101.9-60.0 for figs and 
94.9-56.2 for wasps; Table 1). 

Overall, our analyses favor an Eurasian origin for both 
Ficus and their pollinators. Indeed, in most Eurasian 
clades, Sino-Himalayan figs and their associated 
pollinators appear sister to the rest of the species (Fig. 2, 
gray rhombus). The overall biogeographical signal was 
similar across the different methods used and showed 
instances of dispersal resulting in southward range 
expansion. The major lineages of figs and pollinators 
split during the Tertiary and it appears that they then 
spread southward from Eurasia (Fig. 4), as reflected 
by the branching order of several clades (Fig. 2, gray 
arrows). The major lineages subsequently diversified 
within the Paleotropics and Neotropics during the 
Miocene. 



Discussion 

Codivers ifica Hon 

Our analyses provide both topological and temporal 
lines of evidence to indicate that figs and fig wasps may 
represent the first significant case of long-term (~75 myr) 
codiversification in an insect-plant association. The 
existence of mutualism per se appears insufficient for 
codiversification, because speciation in other intimate 
and sophisticated insect pollination mutualisms (e.g., 
Yuccas and Yucca moths, Glochidion and Epicephala 
moths) seems to be driven by host shifting and host 
tracking rather than cospeciation (Smith et al. 2008; 



Kawakita and Kato 2009). A plausible explanation for 
the significant pattern of cospeciation in the fig-fig 
wasp mutualism is the unusually strong phenotypic 
coadaptation of key traits, such as the specificity of the 
chemical mediation between partners (Grison-Pige et al. 
2002), the lock-and-key shapes of fig ostioles and wasp 
heads (van Noort and Compton 1996; Kjellberg et al. 
2001). 

Despite a history dominated by codiversification, 
there are also some clear mismatches between fig and 
wasp phylogenies (Fig. 2). Our analyses support some 
ancient host-shifts (e.g., by the pollinators of Eriosycea, 
Conosycea, and F. carica), implying that coadapted 
pollinators are sometimes replaced by other wasp 
species without collapse of the mutualism. Finally, 
several host shifts occur at shallow nodes, such 
as between Ficus species in the section Americana 
(Supplementary Fig. S10). 

Overall, our tree reconciliation analyses suggest 
that fig-pollinator history includes numerous species 
duplications and host shifts, as well as cospeciation 
events. However, most host shifts are inferred to have 
occurred between relatively closely related fig species, 
consistent with observations of extant wasp species 
occasionally sharing 2 closely related fig species (Molbo 
et al. 2003; Erasmus et al. 2007). If more distant 
host shifts were common, the congruence of fig and 
wasp phylogenies would be eroded rapidly, even if 
cospeciation remained common (Machado et al. 2005; 
Cook and Segar 2010) . Considering the uncertainty of the 
sister to all other fig-pollinating wasps, it should be noted 
that an alternative topology with Tetrapus as sister to all 
other fig-pollinating wasps would mirror the position of 
Ficus section Pharmacosycea as sister to all other figs and 
should therefore increase cophylogenetic signal. 
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Table 1. Comparison of mean age estimates (Ma) for selected nodes in the fig and wasp phylogenies 



Nodes Ficus mean age Agaonidae mean age 

Ma (95% HPD) Ma (95% HPD) 



Crown Ficus /Crown Agaonidae 


74.9(101.9 


-60.0) 


75.1(94.9- 


-56.2) 


Crown Pharmacosycea/Tetrapus a 
Stem Pharmacosycea/Tetrapus 11 


16.2(25.7- 


-8.2) 


15.9(22.0- 


-9.3) 


74.9(101.9 


-60.0) 


62.1(79.0- 


-45.2) 


Crown Sycomorus/Ceratosolen 


49.1(67.4- 


-34.0) 


59.7(75.4- 


-43.3) 


Stem Sycomorus/ Ceratosolm 


59.4(83.9- 


-43.5) 


69.0(87.6- 


-51.5) 


Crown sect. Adenospermae / 'sg. Strepitus 


35.8(51.9- 


■23.1) 


41.6(57.1- 


-26.7) 


Crown sect. Sycocarpus/sg. Rothropus 


39.5(54.3- 


-26.5) 


A"~t A /[TO 't 

43.4(58.2- 


-30.6) 
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Crown F. pumila group / Wtebesia 
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-11.3) 


32.4(45.3- 


-20.4) 


Stem "F. pumila group" /Wiebesia b 


24.4(36.4- 


- 14.6) 


50.5(63.9- 


-38.1) 


Crown Malvanthera / Pleistodontes 


28.5(41.7- 


-17.8) 


24.1(34.4- 


-14.9) 


Stem Malvanthera/ Pleistodontes 


43.4(61.0- 


29.4) 


50.3(62.0- 


-37.1) 


Crown Oreosycea/Dolichoris 


31.0(46.4- 


-17.6) 


41.0(52.8- 


-30.2) 


Stem Oreosycea/Dolichoris 


46.2(64.6- 


-31.0) 


48.2(61.5- 
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Crown Urostigma/Platyscapa 
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-14.8) 
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-32.7) 


Crown Americana / Pegoscapus 11 
Stem Americana / Pegoscapus* 


20.5(29.3- 
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-28.8) 
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28.3(40.3- 
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Note: 95% lower and upper highest posterior distribution inferred by BEAST is reported between parentheses. 

a Crown-group Pegoscapus and Tetrapus were assigned uniform prior distributions with minimum ages of 15 Ma 

and maximum ages of 60 Ma based on Dominican amber fossils. 

b "F. pumila group" refers to the clade including F. pumila, F. oleifolia, and F. deltoidea. 

c Allotriozoon excepted. 



Biological observations and phylogenetic trees show 
that pollinators of figs are clustered into groups that 
are consistently associated to Ficus sections, subsections, 
and even to some species groups of figs. These inter- 
and intra-generic wasp clusters are highly diverged and 
relatively old and groups of wasps rarely experience 
shifts to other groups of figs. Considering only resolved 
nodes of both phylogenies (Fig. 2, white boxes), we 
observed only 4 shift events between fig subgenera 
[Blastophaga psenes, Wiebesia cf callida, 3 pollinators 
of Frutescentiae {Wiebesia pumilae, W. quadrupes and 
W. sp. ex F. oleifolia) and Valisia spp.] and 5 shift 
events between fig sections [Platyscapa bergi and P. 
sp. (ex F. glaberrima) to Conosycea, Ceratosolen vissali 
and Ceratosolen sp. (ex F. semivestita) to Sycocarpus and 
Adenosperma, respectively, and K. subulatae and K. sessilis 
to Palaeomorphe]. This could be explained by: (i) the 
allopatry of many fig and agaonid groups, (ii) the 
differences between habitats of their host figs (e.g., forest 



canopy versus savannah), and (iii) their host specificity 
due to intricate coadaptation of phenotypes, including 
key traits involved in their reproduction. Consequently, 
we hypothesize that these figs and associated wasp 
groups evolve largely independently as closed systems 
(see also Machado et al. 2005; Cook and Segar 2010) and 
rarely exchange genes or pollinator species, a kind of 
higher level of lineage sorting. 

Recent analyses suggest that the stability of 
fig /pollinator associations can be erratic before 
complete lineage sorting has occurred, or before 
ecological /geographical isolation of fig groups 
(Machado et al. 2005; Jackson et al. 2008; Jousselin 
et al. 2008; Renoult et al. 2009; Cornille et al. 2012). 
During that period of time, pollinator duplication, 
extinction, and hosts shifts within local groups of 
related figs sharing similar phenotypic traits may occur 
frequently. However, these events should not disrupt 
long-term phylogenetic correlations, if lineages sort 
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Figure 4. Hypothetical biogeographical scenario of mutualism diversification. Scenarios are presented on 4 different maps for clarity. Solid 
arrows: figs; dashed arrows: wasps. Nodes with BP < 70% or PP < 0.95 are collapsed. 



over evolutionary time (Cook and Segar 2010). Future 
work should also seek to understand the patterns and 
processes of cospeciation and other processes between 
closely related figs and wasps such as within fig sections. 
Previous studies at this level have focused on sections 
Americana (e.g., Machado et al. 2005; Jackson et al. 2008), 
Galoglychia (e.g., Jousselin et al. 2008), and Sycomorus 
s.l. (Weiblen and Bush 2002) and future studies should 
focus on adding also more dioecious and Eurasian 
clades and to explain why the degree of cospeciation 
appears to vary between ciades. 



Wasp and Fig Systematics 

Our phylogenetic trees provide several new insights 
into the systematics of figs and fig wasps and a sound 
evolutionary framework for future studies in community 
and behavioral ecology. The question of which groups 
of wasps and figs are sister to the rest of agaonids and 
figs, respectively, remains open. Statistical support for 
the deeper nodes of the phylogeny is low and precludes 
us from drawing any definite conclusion. Our additional 
analyses show that there is little support for Tetrapus 
as sister to all other Agaonidae based on molecular 
data (see Appendix S3 in the Supplementary Material 
Online for details). Instead, it appears that Kradibiinae 
(Ceratosolen + Kradibia) or Group 4 (most Wiebesia species 
and pollinators of subsection Frustescentiae) are good 
candidates for the sister taxon to all other agaonids. 
We raise the possibility that an LBA artifact may have 
confounded all previous molecular analyses resulting in 
the inference of Tetrapus as sister to all other agaonids 
(Appendix S3 in the Supplementary Material Online). 

One notable difference between this and previous 
studies concerns taxon sampling, which can be critical 



in phylogenetic analyses (e.g., Philippe et al. 2005). By 
including a larger number of species and using sampling 
that reflects the known diversity of each group, we may 
counteract potential problems with long branches (e.g., 
Bergsten 2005). Increasing taxonomic sampling to break 
up long branches has been applied repeatedly, often 
with the conclusion that earlier studies were misled (e.g., 
Soltis and Soltis 2004; Stefanovie et al. 2004; Leebens- 
Mack et al. 2005; Phillips et al. 2010; Philippe et al. 
2011). 

To further investigate Tetrapus placement, we analyzed 
several morphological characters used in the past 
to assess the relationships between agaonid genera 
(Ramirez 1978; Wiebes 1982b; see Appendix S3 in the 
Supplementary Material Online and Supplementary 
Fig. S14 for details). We show that there is no evidence 
from these morphological characters to support Tetrapus 
as sister to all other agaonids. On the contrary, several 
independent morphological characters support Tetrapus 
as nested within the family and closely related to 
Dolichoris (Appendix S3 in the Supplementary Material 
Online). 

Further studies using more genes and increased 
taxonomic sampling of both ingroups and outgroups of 
figs and wasps are still needed and should contribute to 
resolving the higher taxonomic group relationships with 
more confidence and determine the earliest divergence 
among the figs and their pollinating wasps. 



Pollination 

A number of authors (including some coauthors of 
this study) have previously proposed passive pollination 
as the ancestral mode for the mutualism, followed by a 
single shift to active pollination and several independent 
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reversions to passive (Machado et al. 2001; Jousselin 
et al. 2003; Herre et al. 2008; Jander and Herre 2010). 
This hypothesis has intutive appeal as most other insects 
that pollinate do so passively, but it based primarily 
on the fact that Tetrapus wasps are passive pollinators 
and appeared as the sister of all other pollinators in 
previous phylogenetic analyses. Similarly their host 
figs (Pharmacosycea) appeared as sister to all other figs. 
However, our new phylogenetic trees support a different 
phylogenetic position for Tetrapus. 

Consequently, the issue of ancestral pollination mode 
must be revisited. The phylogenetic tree itself is in 
question, but we also highlight a key issue about 
the interpretation of a given phylogeny. The previous 
conclusion that passive pollination is ancestral relies 
on the assumption that "basal" branches of the trees 
are more informative about ancestral character states. 
However, there is no reason to assume that traits found in 
Tetrapus/Pharmacosycea are "more primitive" or represent 
traits of the common ancestors of both sister groups 
(Krell and Cranston 2004; Crisp and Cook 2005; Lamm 
and Redelings 2009). At the present time, the ancestral 
pollination mode should be considered as equivocal and 
our analyses imply that it remains so. 

Of our 4 reconstructions, 3 find the ancestral state 
equivocal, whereas one (ML on fig phylogeny) favors 
active pollination. This indicates that further studies 
are needed to infer the ancestral pollination mode 
with more confidence. Importantly, these results were 
established on fully bifurcating trees, but in reality the 
backbones of both trees are not strongly supported 
and may change in future studies. Recent advances 
in our understanding of the morphological evolution 
of Moraceae, and in particular of an expanded tribe 
Castilleae, the figs closest relatives, may also shed new 
light on the evolution of the mutualism (Clement and 
Weiblen2009). 



Cobiogeography 

Molecular divergence time estimates point to a 
Cretaceous origin for the mutualism, but differ with 
respect to biogeographic scenarios (Supplementary 
Tables S3 and S5). Previous biogeographic analyses of 
fig- wasps have argued in favor of Gondwanan vicariance 
(Machado et al. 2001). However, a previous study by 
Lopez-Vaamonde et al. (2009) reconstructed ancestral 
areas of fig-pollinating wasps using a phylogenetic tree 
with Tetrapus as sister to the remainder of the fig- 
pollinating wasps. These authors concluded that the 
MRCA of all extant fig-pollinating wasps was most likely 
Asian, although a southern Gondwanan origin could 
not be rejected. A Laurasian origin with subsequent 
dispersal has been proposed for figs and their nearest 
relatives (Zerega et al. 2005). 

Our analyses indicate that the fig-wasp mutualism 
was already in existence ~75 Ma in Eurasia and our 
independently derived mean date estimates for figs 
(75.1 ± 19.4 Ma) and wasps (74.9 ± 21.0 Ma) crown 



groups are remarkably similar, although the size of the 
confidence intervals introduces a degree of uncertainty. 
Despite differences in sampling and dating algorithms, 
the dates obtained correspond well with most previous 
estimates (Supplementary Tables S3 and S5). 

In addition, the hypothesis of an Eurasian origin of 
the mutualism is supported by several other lines of 
evidence: (i) the presence in Asia of 70% of the major 
Ficus clades; (ii) the early divergence of Sino-Himalayan 
fig and wasp species in most Eurasian clades (e.g., F. 
henryi, F. sarmentosa, F. tikoua, F. nervosa); (iii) the fact 
that pollinators of the subsection Frustescentiae are found 
only in Continental Asia (Fig. 2); (iv) the age estimates 
for Moraceae in general and Ficeae (Dorstenieae and 
Castilleae) in particular (Zerega et al. 2005); and (v) 
the fact that the oldest fig and wasp fossils are known 
only from the Northern Hemisphere (Collinson 1989; 
Compton et al. 2010) (see our review of the literature on 
fig fossils in Appendix S4 in the Supplementary Material 
Online). Finally, Burnham and Graham (1999), analyzing 
the origin of the tropical component in northern Latin 
American vegetation also suggest that Ficus arrived 
from the north. Accordingly, current data support the 
conclusion that the mutualism probably originated 
in the tropical forests of Eurasia (Otto-Bliesner and 
Upchurch 1997). 

Pharmacosycea and Tetrapus divergence is dated to 74.9- 
62.1 Ma (mean stem figs-mean stem age wasp; Table 1), 
before South America split from Antarctica. However, 
rather than explaining the South American colonization 
of Pharmacosycea /Tetrapus by trans-antarctic routes, we 
propose that both lineages might have reached the 
New World across North Atlantic land bridges (Tiffney 
1985), dispersing through the evergreen woodland and 
tropical forest belts of Eurasia (Fig. 4). South America 
may have been colonized later via "stepping-stone" 
volcanic islands. Indeed, most Pharmacosycea species 
inhabit the Northern Andes and there are none in Chile 
and Patagonia, which have vegetation similar to late 
Cretaceous Antarctica (Poole et al. 2003). Because figs 
are also absent from the exceptionally good fossil flora 
of Patagonia (Wilf et al. 2003), trans-Antarctic dispersal 
seems unlikely but cannot be completely ruled out. 
Although the Laguna del Hunco flora is not considered a 
tropical flora, this Patagonian flora hosts one Papuacedrus 
species very closely related to extent tropical Papuan 
species (Wilf et al. 2009). In Papua New Guinea and in 
Papua Barat (Indonesia), this conifer is found in the same 
mountainous habitats as Malvanthera fig tree species at 
altitudes > 2000m (for example in the Arfak mountains, 
Kebar Valley, Bulolo-Wau, Mt Kerewa). Accordingly, 
despite not being considered a tropical flora, the Laguna 
del Hunco flora could have also have included Ficus and 
the fact that Ficus appears absent from this exceptionally 
good flora, supports the later arrival of Ficus from 
Eurasia. 

Based on our biogeographic analyses, the major 
lineages of figs and pollinators split during the Tertiary 
and spread southward from Eurasia (Fig. 4), possibly 
in response to the cooling climate (Davis et al. 2002). 
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Subsequent diversification occurred within continents 
during the warmer Miocene epoch (Zachos et al. 2001). 

The general scenario of fig-wasp codiversification 
is illustrated by the charismatic hemi-epiphytic or 
"strangler" figs (the subgenus Urostigma clade, Fig. 4), 
which evolved ~52-50.3 Ma, during a period of global 
warming. A first clade, probably living west of the Turgai 
straits (Akhmetiev and Beniamovski 2009), dispersed 
southward into Africa to form section Galoglychia and 
into South America to form section Americana, some 
32.3-38.2 Ma. Another clade, probably occurring in 
east Eurasia, spread to India and Sundaland to form 
section Conosycea and to Australasia to form section 
Malvanthera, ~50.3-43.4 Ma. This latter dispersal was 
probably via stepping-stones through the Ninety East 
Ridge (Carpenter et al. 2010), because direct dispersal 
from Sundaland to Australia was impossible before 
25 Ma (Hall 2002). Today, each tropical continent has its 
own major endemic radiation of strangler figs, stemming 
from these ancient dispersal processes. Interestingly, 
pollinator biogeography shows a few discrepancies 
with this scenario for fig dispersal. Indeed, the genus 
Pleistodontes pollinating Malvanthera figs is sister to all 
other Urostigma pollinators. Therefore, we propose that 
Conosycea was colonized by a host shift of an ancestral 
Galoglychia/ Americana pollinator in southern Eurasia 
before spreading to southern Sunda. 



Supplementary Material 

Supplementary material, including data files 
and /or online-only appendices, can be found in 
the Dryad data repository at http://datadryad.org, doi: 
10.5061 /dryad.hr620. Matrices are also available 
in TreeBASE (No. TB2:S13315) at http:/ /www. 
treebase.org/. 
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Conclusion 

Based on multiple lines of evidence (fossils, 
Moraceae history, branching pattern, and ancestral 
area reconstructions), we infer an Eurasian origin 
for the fig /pollinator mutualism. We show that the 
mutualism arose ~75 Ma, confirming previous estimates 
(Supplementary Tables S3 and S5). Because that time, 
the insects and plants have diversified together 
leaving a strong long-term signal of phylogenetic 
congruence, confirming previous studies based on 
smaller data sets (Supplementary Tables S3 and S5). 
This is not due to strict cospeciation alone, but reflects 
a history with large amounts of cospeciation and 
insufficient host shifts to alter the marked phylogenetic 
matching. This is the only known example of long- 
term insect /plant codiversification and we are not 
aware of other candidates for such a pattern. Other 
insect/plant pollination mutualisms do not appear to 
be characterized by phylogenetic congruence, and we 
propose that strong codiversification of figs and their 
pollinators is driven by their unusually high level of 
phenotypic trait matching. Figs and their pollinators 
have spread across the globe to occupy all tropical 
continents, where they play important ecological roles 
in forests and savannahs. Their numerous interactions 
with other species, such as vertebrate frugivores, mean 
that the evolution of entire tropical ecosystems has been 
influenced strongly by this unique strong pattern of 
codiversification between fig trees and their pollinating 
insects. 
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